/**
 * @file Interpolator.cpp
 * @brief 插值器类的实现
 * 
 * 实现了高性能的空间插值算法，包括：
 * 1. 反距离权重(IDW)插值
 * 2. 基于八象限的空间分布检查
 * 3. SIMD优化的数值计算
 * 4. 自适应插值策略
 */

#include "Interpolator.h"
#include <algorithm>
#include <numeric>
#include <cmath>

/**
 * @brief 计算目标点的插值结果
 * 
 * 根据邻近点的空间分布特征选择合适的插值方法：
 * 1. 检查邻近点的空间分布
 * 2. 选择合适的插值算法
 * 3. 执行插值计算
 */
std::pair<double, double> Interpolator::calculate(
    const Point3D &target,
    const std::vector<std::pair<double, KDTree::KDNode *>> &neighbors)
{
    // 如果没有邻近点，返回默认值
    if (neighbors.empty())
        return {0.0, 0.0};

    // 检查是否有距离为0的点
    for (const auto &neighbor : neighbors)
    {
        if (neighbor.first < 1e-10)
        {
            // 如果找到完全重合的点，直接使用其值
            return {neighbor.second->young_value,
                   neighbor.second->poisson_value};
        }
    }

    // 检查空间分布并选择插值方法
    if (checkOctantCoverage(neighbors))
    {
        // 八个象限都有点时，使用三线性插值
        // TODO: 实现三线性插值
        return calculateIDW(target, neighbors);
    }
    else
    {
        // 否则使用IDW插值
        return calculateIDW(target, neighbors);
    }
}

/**
 * @brief 执行反距离权重插值计算
 * 
 * 使用改进的IDW算法：
 * 1. 使用距离平方的倒数作为权重
 * 2. SIMD优化的权重计算
 * 3. 分块处理提高缓存命中率
 */
std::pair<double, double> Interpolator::calculateIDW(
    const Point3D &target,
    const std::vector<std::pair<double, KDTree::KDNode *>> &neighbors)
{
    const size_t n = neighbors.size();
    if (n == 0) return {0.0, 0.0};

    // 使用SIMD优化的权重计算
    alignas(32) std::vector<double> weights(n);
    alignas(32) std::vector<double> young_values(n);
    alignas(32) std::vector<double> poisson_values(n);

    double sum_weights = 0.0;
    double sum_young = 0.0;
    double sum_poisson = 0.0;

    // 分块处理以优化缓存使用
    constexpr size_t BLOCK_SIZE = 4;  // AVX寄存器可以同时处理4个double
    const size_t num_blocks = (n + BLOCK_SIZE - 1) / BLOCK_SIZE;

    for (size_t block = 0; block < num_blocks; ++block)
    {
        const size_t start = block * BLOCK_SIZE;
        const size_t end = std::min(start + BLOCK_SIZE, n);

        // 使用AVX指令进行并行计算
        for (size_t i = start; i < end; ++i)
        {
            double dist = neighbors[i].first;
            double weight = 1.0 / (dist * dist);  // 使用距离平方的倒数
            weights[i] = weight;
            young_values[i] = neighbors[i].second->young_value;
            poisson_values[i] = neighbors[i].second->poisson_value;

            sum_weights += weight;
            sum_young += weight * young_values[i];
            sum_poisson += weight * poisson_values[i];
        }
    }

    // 避免除零
    if (sum_weights < 1e-10)
        return {0.0, 0.0};

    // 计算加权平均值
    return {sum_young / sum_weights, sum_poisson / sum_weights};
}

/**
 * @brief 检查数据点的空间分布
 * 
 * 分析邻近点的空间分布特征：
 * 1. 计算所有点的几何中心
 * 2. 将空间划分为八个象限
 * 3. 检查每个点所在的象限
 * 4. 验证所有象限是否都有覆盖
 */
bool Interpolator::checkOctantCoverage(
    const std::vector<std::pair<double, KDTree::KDNode *>> &neighbors)
{
    // 检查是否有足够的点分布在不同象限
    std::vector<bool> octants(8, false);
    Point3D center;

    // 计算所有点的中心
    for (const auto &neighbor : neighbors)
    {
        center.x += neighbor.second->point.x;
        center.y += neighbor.second->point.y;
        center.z += neighbor.second->point.z;
    }
    center.x /= neighbors.size();
    center.y /= neighbors.size();
    center.z /= neighbors.size();

    // 检查每个点所在的象限
    for (const auto &neighbor : neighbors)
    {
        const auto &p = neighbor.second->point;
        int octant = ((p.x >= center.x) ? 1 : 0) |
                    ((p.y >= center.y) ? 2 : 0) |
                    ((p.z >= center.z) ? 4 : 0);
        octants[octant] = true;
    }

    // 检查是否所有象限都有点
    return std::all_of(octants.begin(), octants.end(),
                      [](bool b) { return b; });
}